# set working directory to location of the R file
#setwd("replication/scripts")

require(devtools)
# uncomment if not installed in previous file:
#devtools::install_github("ChandlerLutz/starpolishr")
library(starpolishr)

require(pacman)
pacman::p_load(tidyverse, lfe, dplyr, ggplot2, ggpubr, stargazer,
               purrr, cowplot, starpolishr, here)

# All models run using PanelMatch version 1.0.0

require(devtools)
# uncomment if not installed in previous file:
#install_version("PanelMatch", version = "1.0.0", repos = "http://cran.us.r-project.org")

library(PanelMatch)

# combined panel
load("../Data/combined_panel.RData")

load("../Data/staff_party.RData")

staff_firms_party_panel <- select(staff_firms_party_panel, reg_lobby_core_id, period, any_dem, any_majority, personal_staff)
combined_panel <- left_join(combined_panel, staff_firms_party_panel)

year_periods <- data.frame(period = seq(1,42), year = rep(seq(2000, 2020), each=2))

combined_panel <- combined_panel %>% 
  left_join(year_periods) %>% 
  group_by(reg_lobby_core_id, year) %>% 
  fill(any_majority, .direction = "updown") %>% 
  fill(personal_staff, .direction = "updown")


combined_panel <- combined_panel %>% 
  ungroup %>% 
  mutate_at(vars(any_dem, any_majority), ~replace_na(., 0)) %>% 
  mutate(any_minority = ifelse(any_majority == 0 & personal_staff == 1, 1, 0) %>% replace_na(0)) 





# bureaucrat and staff separate panels
load("../Data/firm_experience_panels.RData")

###
# bureaucrats
prev_exp_firm <- as.data.frame(prev_exp_firm)
prev_exp_firm$period <- as.integer(prev_exp_firm$period)
prev_exp_firm$log_rev <- log(prev_exp_firm$revenue)

# check staff and bureau covariates are the same
prev_exp_firm$log_contracts <- log(prev_exp_firm$total_filings)
prev_exp_firm$log_lobby <- log(prev_exp_firm$total_lobbyists)

# create treatment indicator
prev_exp_firm$prev_exp2 <- ifelse(is.na(prev_exp_firm$prev_exp)==T,0,prev_exp_firm$prev_exp)

###
# staff
staff_exp_firm <- as.data.frame(staff_exp_firm)
staff_exp_firm$period <- as.integer(staff_exp_firm$period)
staff_exp_firm$log_rev <- log(staff_exp_firm$revenue)

# check staff and bureau covariates are the same
staff_exp_firm$log_contracts <- log(staff_exp_firm$total_filings)
staff_exp_firm$log_lobby <- log(staff_exp_firm$total_lobbyists)

# create treatment indicator
staff_exp_firm$staff_exp2 <- ifelse(is.na(staff_exp_firm$staff_exp)==T,0,staff_exp_firm$staff_exp)



#############
# convenience functions

# function for extracting id of control firms

extract_control <- function(cont, pm_obj){
  
  crnt_control <- as.numeric(as.character(pm_obj$att[[cont]]))
  return(crnt_control)
}

# wrapper function for PanelMatch

pm_wrap <- function(lags, dataset, treat, refine = F,
                    DV){
  
  if(refine == F){
    comb <- PanelMatch(lag = lags, 
                       time.id = "period", unit.id = "reg_lobby_core_id", 
                       treatment = treat, 
                       refinement.method = "none", 
                       data = dataset, 
                       match.missing = T, 
                       qoi = "att" ,
                       outcome.var = DV,
                       lead = 0:4, 
                       forbid.treatment.reversal = T)
    
    dist_control <- comb$att
    comb_obs <- sum(print(comb$att)[,3])
    
    if(nchar(comb_obs) == 6){
      comb_obs <- paste(substr(comb_obs, 1, 3), ",", 
                        substr(comb_obs, 4, nchar(comb_obs)), sep = "")}else{
                          comb_obs <- paste(substr(comb_obs, 1, 2), ",", 
                                            substr(comb_obs, 3, nchar(comb_obs)), sep = "")
                        }
    
    unique_treat <- length(unique(print(comb$att)[,1]))
    
    map_length <- 1:length(print(comb$att)[,1])
    unique_control <- purrr::map(map_length,
                                 ~ extract_control(cont = ., pm_obj = comb))
    
    unique_control <- do.call("rbind", unique_control)
    unique_control <- length(unique(unique_control))
    
    unique_control <- paste(substr(unique_control, 1, 2), ",", 
                            substr(unique_control, 3, nchar(unique_control)), sep = "")
    
  }else{comb <- PanelMatch(lag = lags, 
                           time.id = "period", unit.id = "reg_lobby_core_id", 
                           treatment = treat, 
                           refinement.method = "mahalanobis", #alternative none
                           data = dataset, 
                           match.missing = T, 
                           covs.formula = ~ I(lag(log_contracts,1:4)) + I(lag(log_lobby,1:4)) + I(lag(log_rev,1:4)), 
                           qoi = "att" ,
                           outcome.var = DV,
                           lead = 0:4, 
                           forbid.treatment.reversal = T)
  dist_control <- comb$att
  comb_obs <- sum(print(comb$att)[,3])
  
  if(nchar(comb_obs) == 6){
    comb_obs <- paste(substr(comb_obs, 1, 3), ",", 
                      substr(comb_obs, 4, nchar(comb_obs)), sep = "")}else{
                        comb_obs <- paste(substr(comb_obs, 1, 2), ",", 
                                          substr(comb_obs, 3, nchar(comb_obs)), sep = "")
                      }
  
  unique_treat <- length(unique(print(comb$att)[,1]))
  
  map_length <- 1:length(print(comb$att)[,1])
  unique_control <- purrr::map(map_length,
                               ~ extract_control(cont = ., pm_obj = comb))
  
  unique_control <- do.call("rbind", unique_control)
  unique_control <- length(unique(unique_control))
  
  unique_control <- paste(substr(unique_control, 1, 2), ",", 
                          substr(unique_control, 3, nchar(unique_control)), sep = "")
  
  
  }
  
  comb_results <- PanelEstimate(sets = comb, 
                                data = dataset, 
                                confidence.level = .9)
  
  pm_res <- as.data.frame(summary(comb_results)[[1]][1:5, 1:2])
  
  
  return(list(res = pm_res, 
              n_control = comb_obs, 
              n_treat = unique_treat,
              n_unique_control = unique_control,
              dist_control = dist_control))
  
}


combined_panel <- as.data.frame(combined_panel)
combined_panel$period <- as.integer(combined_panel$period)
combined_panel$log_rev <- log(combined_panel$revenue)

combined_panel$log_contracts <- log(combined_panel$total_bur_filings)
combined_panel$log_lobby <- log(combined_panel$total_bur_lobbyists)

combined_panel$treat <- ifelse(combined_panel$prev_exp == 1 | combined_panel$staff_exp == 1, 1, 0)
combined_panel$treat <- ifelse(is.na(combined_panel$treat)==T,0,combined_panel$treat)

comb1 <- pm_wrap(lags = 6, dataset = combined_panel,
                treat = "any_majority", refine = T, DV = "log_rev")

comb2 <- pm_wrap(lags = 6, dataset = combined_panel,
                treat = "any_minority", refine = T, DV = "log_rev")



mod <- felm(log_rev ~ treat
            |reg_lobby_core_id + period | 0|reg_lobby_core_id,
            data = combined_panel) 


# Panel A
pan1 <- stargazer(mod, mod, mod, mod, mod, 
                  coef =  comb1$res[,1], se = comb1$res[,2],
                  covariate.labels = c("Lobbyist as Staffer (in Majority)"),
                  column.labels = c("t=0", "t+1", "t+2", "t+3", "t+4"), no.space = T,
                  omit.stat = c("f", "adj.rsq", "rsq", "ser", "n"),
                  dep.var.labels = "ln Revenue", 
                  title = c("Lobbying Firm Revenue and Majority versus Minority Connections"), 
                  label = c("majority_minority"),
                  add.lines = list(c("Treatment Events (Majority)", comb1$n_treat, comb1$n_treat, comb1$n_treat,
                                     comb1$n_treat, comb1$n_treat),
                                   c("Treatment Events (Minority)", comb2$n_treat, comb2$n_treat, comb2$n_treat,
                                     comb2$n_treat, comb2$n_treat)),
                  header = F)

#Panel B
pan2 <- stargazer(mod, mod, mod, mod, mod, 
                  coef =  comb2$res[,1], se = comb2$res[,2],
                  covariate.labels = c("Lobbyist as Staffer (in Minority)"),
                  column.labels = c("t=0", "t+1", "t+2", "t+3", "t+4"), no.space = T,
                  omit.stat = c("f", "adj.rsq", "rsq", "ser", "n"),
                  dep.var.labels = "ln Revenue", 
                  title = c("Lobbying Firm Revenue and Majority versus Minority Connections"),
                  #label = c("combined_adj"),
                  # add.lines = list(c("Treatment Events", comb2$n_treat, comb2$n_treat, comb2$n_treat,
                  #                    comb2$n_treat, comb2$n_treat),
                  #                  c("Control Firms", comb2$n_control, comb2$n_control,
                  #                    comb2$n_control, comb2$n_control, comb2$n_control),
                  #                  c("Unique Controls", comb2$n_unique_control, comb2$n_unique_control,
                  #                    comb2$n_unique_control, comb2$n_unique_control, comb2$n_unique_control)),
                  header = F)

# put combined table together
BothPanels <- star_panel(pan1, pan2, panel.names = c("Majority Staffer Connection", "Minority Staffer Connection"))

# add a note
BothPanels <- star_notes_tex(BothPanels, note.type = "threeparttable",
                             note = "Estimates are the averages of all possible 2x2 differences-in-differences where control groups include firms with identical treatment histories prior to each respective 
                             treatment event. Future-treated firms and treatment-switchers are never included. Panel A includes lobbyists as staffers in personal offices in the majority party; Panel B includes lobbyists as staffers in personal offices of the minority party.
                             All models' control groups are made by matching firms using the Mahalanobis 
                             distance calculated on logged firm revenue, logged number of contracts and logged number of lobbyists, all calculated in the pre-treatment windows. Weighted firm-blocked
                             bootstrapped standard errors in parentheses. 1,000 trials used. *p$<$0.1; **p$<$0.05; ***p$<$0.01")

# export
star_tex_write(BothPanels, 
               file = "../tables/TableD11.tex")



